: A Version of Simpson's Rule for Multiple 



bjol Integrals 

^3 ■ 



< 
00 



Alan Horwitz 
February 8, 2008 



< 

: 1 Introduction 

Let M(/) denote the Midpoint Rule and T(/) the Trapezoidal Rule for esti- 
mating JI^ f[x)dx. Then M{f) and T(/) are each exact for f{x) = 1 and x. 
A more accurate rule can be obtained by taking XM{f ) + (1 — A)T(/), where 
ly-^ I A = |. Of course this rule is known as Simpson's Rule, and is exact for all 

^ ■ polynomials of degree < 3. The purpose of this paper is to extend this idea 

QQ ■ to multiple integrals over certain polygonal regions D„ in i?". First, consider 

O . the case n = 2. So suppose that Z) is a polygonal region in R^. The midpoint 

rule in one dimension is M{f) = (h — a)f (^^^^ • A natural extension of this 
rule is 

M(/) = A{D)f\P) 



0^ 



where A{D) = area of D, and P = centroid of D. 

The trapezoidal rule in one dimension is T(/) = (h ~ a) |^/MtiM^ . A 
^ \ natural extension of that rule is 

^' T{f) = AiD)-f:m) 

where the Pk are the vertices of 

In general, let D„ be some polygonal region in _R", let Pq, ...,Pm denote 
the vertices of and let Pm+i = center of mass of Dn- Define the linear 
functionals 



° Key Words: Cubature Rule, Simpson's Rule, polygonal region, Grobner Basis, exact 
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M(/) = Vol(D„) f{Pm+i), nf) = Vol(D„)(^E7=o/(^i)), and for 
fixed A, < A < 1, 

L, = XM{f) + (1 - A)T(/) 

The idea is to choose A so that Lx is a good cubature rule(CR). Our objective 
here is not to provide optimal cubature rules, but to generalize the ideas 
from Simpson's Rule in one variable to several variables. In some cases we 
have reproduced known CRs using a different approach^. In other cases our 
CRs appear to be new. In addition, our approach suggest a general method 
for deriving CRs for polygonal domains. 

For the n cube [0, 1] x ■ ■ ■ x [0, 1], A = |, and the CR L{f) is exact for 
all polynomials of degree < 3 (see CRS). This, of course, matches the degree 
of exactness of Simpson's Rule in one variable. For the n simplex, A = 
and the CR L{f) is exact for all polynomials of degree < 2(see CRl). 

For regular polygons in the plane with m sides, our approach fails for 
m > 4. For example, if D = a. regular hexagon, then no choice of A makes 
Lx exact for all polynomials of degree < 2. Even for m = 4, if the polygon 
is not regular, then our approach fails as well. For example, let D be 
the trapezoid with vertices {(0,0), (1,0), (0, 1), (1,2)}. Again, no choice of 
A makes Lx exact for all polynomials of degree < 2. One can do better, 
however, by using points Qj G d{Dn), other than the vertices of Dn, to 
generate T{f). In some cases this leads to better formulas. For example, for 
the trapezoid D, this leads to a formula which is exact for polynomials of 
degree < 2 (see CR 5). Choosing points different from the vertices leads to 
a system of polynomial equations. We use Grobner Basis methods(see P), 
along with Maple, to solve such systems when possible. For the n Simplex 
Tn, we try using the center of mass of the faces of T„ to generate T(/)(see 
CR2). A similar idea works for the n cube. 

Instead of using the weighted combination, XM{f) + (1 — A)T(/), another 
way to derive Simpson's Rule in one variable is to integrate the quadratic 
interpolant to / at a, and b. For the n Simplex T„, our generalization 
Lx can be obtained by integrating p{x) over T„, where p{x) is the unique 
interpolant to f{x) at the vertices and center of mass of T^. Indeed, M(/) = 
Jrp^ T{x)dV, where T{x) = tangent hyperplane to /(x) at the center of mass 
of Tn, and T(/) = /^^ L{x)dV, where L{x) = bilinear interpolant to f{x) at 
the vertices of T„. This, of course, is the generalization from one variable, 

"'^See and for a thorough treatment of CRs. 
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where {b - a)f (^) = J^T{x)dx, T = tangent line to f(x) at and the 
trapezoidal rule can be obtained by integrating the linear interpolant to / 
at the endpoints. For regions -D„ in general, however, Lx,T, and M do not 
always arise in this fashion. Indeed, the number of basis functions does not 
always match the number of nodes, and the interpolant may not be unique. 
This is what occurs for the n cube. 

We also give an analgous formula, using four knots, for the unit disc(see 
CR 6). 

Finally we note that for most of our CRs, all of the weights are positive, 
all of the knots he inside the region D„, and all but one of the knots hes on 
the boundary of the region. 

2 n Simplex 

Let Pq, ...,Pn denote the n + 1 vertices of the n simplex T„ C i?". Hence 
Pj = (0,0, ...,1,0,0) for k >1, and Pq = (0, 0, 0), x = {xi, ...,Xn),dV = 
standard Lebesgue measure on T„. We hst the following useful facts: 

Vol(r„) = ^, P„+i = Center of Mass of r„ = (l/(n + 1), l/(n + 1)) 
nl 

l^x,dV=^^,k^l,...,n 
r 2 

/ xldV = —,k — l,...,n 

JTr. ^ (n + 2)!' 

Define the following linear functionals: 

M(/) = Vol(r„) /(P„+i) = -f{l/{n + 1), l/(n + 1)) 

77/. 

-in -in 

T(/) ^ voi(r„,(— g/w)) ^ ^^^J(P>) 

Hf) = / f(x)dV 
Finally, for fixed A, < A < 1, define 
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= AM(/) + (1 - A)r(/) 

First, if fix) = Xk, then /(/) = M(/) = T{f) = for any A. Now 

let fix) = xjxk, J ^ k. Then /(/) = jj^,,M{f) = ^-^^X__,t(/) = 0. 
Hence = /(/) ^ A ((TttTiW) = (^IW ^ ^ = Let L = L^.with 

this A. Now, if fix) = xl then L(/) = ^^(^^ + ^^(^^ = = 
/(/). We summarize 

72+1 1 

CRl: Lif) = - — — — /(P„+i)+- — — — - V fiPj) is exact for all polynomials of degree < 2. 
in + 2)n\ ("- + 2)!~^ 

L is not exact for all polynomials of degree < 3. For example, if fix) = 
XjXkXiJ ^ k,k ^ l,j ^ /, then /(/) = while L(/) = (^^^(^^2)' ^ 

/(/)• 

2.0.1 Connection with Interpolation 

Another way to derive Simpson's Rule in one variable is by using quadratic 
interpolation at a, b. The cubature rule Lif) above can also be obtained 
by integrating a certain second degree interpolant to / over T„. Let 

n 

piXi, ...,Xn) = An+lXiX2 + AkXk + Aq 

k=l 

We wish to choose the Aj so that 

piPk) = fiPk), k = 0,...,n+l 



It follows easily from (|]0]T|) that Aq = /(Po),^fc = fiPk),k = 
and 

An+i = fiPn+i) -in + l) Yrk=o fiPk). Hence pix)dV = A^+^j^ + 
j^, ELi + h.^o = j^. {{n + 1) -in + l) ELo f{Pk)) 

+ {EU fiPk) - nfiPo))+^fiPo) = ^m+i) + j;^ ELo fiPk) = 
Lif). 

Remark 1 Note that one could use any XjXk,j ^ k in constructing the 
interpolant p, and still obtain JT„pix)dV = Lif). 
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Remcirk 2 This, of course, generalizes the fact that Simpson's Rule in one 

variable can be obtained by integrating the quadratic interpolant to f at 
a,b, In more than one variable, however, our quadratic interpolant only 
uses one second degree basis function. Using the basis functions {1, xi, ...,Xn, X1X2}, 
one obtains exactness for those functions as well as for " "^^""^ additional 
quadratic terms. 

Remcirk 3 One can also obtain the rules M{f) and T{f) using interpola- 
tion. 

M(f) — Jj,^ T{x)dV, where T is the tangent hyperplane to / at Pn+i- 
T{f) = It„ Q{x)dV, where q{xi, Xn) = ELi ^kXk + Bq, the Bk chosen 

so that 

g(Pfc) =/(Pfe), k = 0,...,n. 

Again, this generahzes the fact that the midpoint rule in one variable 
can be obtained by integrating the tangent line at the midpoint, while the 
trapezoidal rule can be obtained by integrating the linear interpolant to / at 
the endpoints. 

2.1 Other Boundary Points 

It is interesting to examine what happens if we use other points on d{Sn) to 
generate T(/). We first examine n = 2 in some detail, and then generalize. 

2.1.1 n = 2 

Let T2 = the triangle in with vertices {(0, 0), (1, 0), (0, 1)}. Use the points 
(a, 0), (0, 6), (c, 1 - c), < a < 1, < 6 < 1, < c < 1 from each side of T2 to 
define 

T(/) = Area(r2)^(/(a, 0)+/(0, 6)+/(c, 1-c)) = l(/(a, 0)+/(0, 6)+/(c, 1-c)) 
As earlier (with (1/3, 1/3) = Center of mass of T), 

M(/) = Area(r2)/(l/3, 1/3) = 1/(1/3, 1/3) 
and, for fixed A, < A < 1, 

= AM(/) + (1 - A)T(/) 
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. L^{x) = iA + |(l-A)(a + c) 
. L,(y) = iA + |(l-A)(6+l-c) 
. L,{x') = ±X + l(l-X)(a' + c') 
. L,(|/^) = ^A + |(l-A)(62 + (l-c 
. L,(xy) = ^A + |(l-A)c(l-c) 

We want to choose a, b, c, and A so that Lx is exact for the functions 
X, ?/^,and xy. Setting Lx{f) = I{f) yields the system of equations 
iA + i (1 - A) (a + c) = |, |A + 1(1 - A) (6 + 1 - c) = |, 
3^A + |(1-A) (a2 + c2) = ^, ^A + |(l-A) (fe^ + (1 - ) = ^, 
^A + |(l-A)c(l-c) = ^ 

We found the following Grobner basis using Maple: 

{-I2c^ + 12Ac^ + 4A - 3 + 12c - 12Ac, -l + a + c,b-c} 

Setting each polynomial from the Grobner basis to yields precisely the same 

solutions as the original system. 

Hence 6 = c, a = 1 - c, and A = ^g^^^, 1 - A = ijg^ 

This gives Lxif) = Ai/(l/3, l/3) + (l-A)i(/(l-c, 0)+/(0, c)+/(c, 1-c)) 

Now we choose c so that Lx is exact for quadratics. Since Lx{x^) = 

^9c!^i±3c^^gg^ and x'dA = ^, Lxix') = ± ^ c = or c = c = Q 

uses the vertices of T2 to generate T(/). Hence we use c = | ^ A = 0. 
With this choice of c and A, it follows that Lxiy"^) = J-p^ y^dA and Lx{xy) = 
Jt^ xydA. 

We summarize 

Lif) = ;^(/(^,0)+/(0, h+f{^, h) is exact for all polynomials of degree < 2. 
o 2 III 



Remark 4 The CR above is given in the well known book of Stroud (see 
as one of a group of formtdas for the triangle T2. 

Since L(x^) = ^ 7^ /tj x^dA = ^, L is not exact for all cubics. 
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Connection with Interpolation It is not hard to show that there is a 
unique interpolant, p(x, y) = Aixy + A2X + A^y + A4, to any given /(x, y) at 
(1/2, 0), (0, 1/2), (1/2, 1/2), (1/3, 1/3). It also follows easily that /^^ p(x, y)dA = 
|(/(i, 0) + /(O, i) + /(I, I)), which is L{f) given above. 

2.1.2 General n 

For T2 above we used the midpoints of the edges for r(/), so it is natu- 
ral to try using the center of mass of the faces of T„, given by Qk — 
^, 0, ^)(all coordinates ^ except the Mh coordinate = 0), A; = 
1, n, Qn+i = (^, ^, ^). This yields 

1 n+1 

n/) = (;^i:/(ft) 

As earlier, 

M(/) = ^/(l/(n + l),...,l/(n+l)) 

= AM(/) + (1 - A)r(/) 

A simple computation shows that Lx is exact for Xk for any A, and is exact 
for x^x,,j ^ k, if (-1_)A+ (ii^)(l - A) = This gives A = - (n - 2) 
With this A, we have 

CR2: L(/) = - (n - 2) V(^+l))+(;^ E /(Q^) 

is exact for all polynomials of degree < 2. Note that one of the weights is 
< if n > 3. 

Remcirk 5 It would be nice to choose points on the faces which yield A > 0, 

and in particular A = 0. We tried this for n = 3, using points on each 
face of the simplex T3. Letting Qi = (oi, 02, 0), = (03, 0, 04), = 
(0, as, Og), Qa = {cir, cis-, 1 — a? — a^), with 

dj > Vj, aj + Cj+i < 1, j = 1, 3, 5, 7 



7 



T(/) = area(T3)lx:/(^?.) = ^E/(^?.) 

M(/) = i/(l/4,l/4,l/4) 

= AM(/) + (1 - X)T(f) 

The equations giving exactness for x, y, z, x^, y^, z'^, xy, xz, yz, after some 
simphfication, look like ai + 03 + oy = 1 , 02 + 05 + ag = 1 , 04 + ag — 07 — Og = 0, 
^ A + (1— A)^(aia2 + ^tO-s) = 
^ A + (1 - A)^(a3a4 + 07(1 - aj - ag)) = 
^ A + (1- A)^(a5a6 + 03(1 - ay - ag)) = j^, 
^ A + (l- X)^^{al + al + af) = ^, 

gg A + (1- A)24(a2 + '^5 + = 60' 

^ A + (1- A)^(a2 + al + {l-ar- ag)^) = ± 

One can solve the first two equations for aj and Og, and then substitute 
into the rest of the equations, yielding seven polynomial equations in seven 
unknowns. Of course two solutions to these equations are 

• ai = 02 = 04 = ae = 0,7 = ag = 0, = = 1, A = |, which uses the 
vertices of T3 for Qj , and 

• = I for all J, , A = — |, which uses the center of mass of each face 
of T3 for Qj. 

We applied Grobner Basis methods to the above equations giving exact- 
ness for X, y, z, x'^, y'^, z'^, xy, xz, yz, along with A = 0. We were then able to 
show that there is no solution. That proves that one cannot make A = 0. 
We were not able to find any positive values for A, though it is possible that 
some exist. 

3 Unit n Cube 

Let C„ = [0, 1]" C -R", and let Pq, Pm-i denote the m vertices of m = 
2"^. Assume that Pq = (0,...,0) and Pm-i = (I,---,!)- The other vertices 
have at least one coordinate which equals 0, and at least one coordinate 
which equals 1. Let x = {xi, ...,Xn),dV = standard Lebesgue measure on 
C„. We hst the following useful facts: 
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vol(C„) = 1, Center of Mass = Pm = (1/2, . . . , 1/2) 
For any non-negative ri, r„, 

/ ^kAl---4'dV= V V ^— 

Jc^ ri + 1 r2 + 1 Tj + 1 

Define the following linear functionals: 

M(/) = Vol(C„) f{Pj = /(1/2, 1/2) 

-1 m— 1 -1 m— 1 

T(/) = voi(c„)(- E /(^.O) = - E /(^.O 

I{f) = / /(S)^^^ 

Finally, for fixed A, < A < 1, 

L, = XM{f) + (1 - A)T(/) 

First, let /(x) = x^- Then M(/) = i, and it is not hard to show that 
T(/) = i2"-i = i. Hence M(/) = T(/) = /(/) ^ L^if) = /(/) for any A. 
Now let f{x) = XjXkij 7^ k. Then M(/) = i. since there are 2""^ ways to get 
a 1 in both the jth and /ct/i coordinates of a vertex of D, T(f) — ^2""^ = |. 
Hence, again, M(/) = T(/) = /(/) for any A ^ La(/) = /(/). Also, if 
fix) = xl then M(/) = 1 and T(/) = 1. Hence L;,(/) = /(/) = i ^ 
A = |. Now let L = L2/3. Finally we consider third degree terms. First, if 
fix) = x,XkXi,j ^k,k^l,j^l, then T(/) = M(/) = /(/) = | ^ L(/) = 
/(/) for any A. Second, if fix) = x]xk,j 7^ A;, then T(/) = i, M(/) = |, and 
/(/) = i, and hence |M(/) + iT(/) = /(/). Finally, if fix) = xl then 

T(/) = |,M(/) = |, and /(/) = i, and again |M(/) + |r(/) = /(/). 
We summarize 

2 1 1 2"-^ 

CR3: L(/) = -/(1/2, 1/2)+ — - J2 f{Pj) is exact for all polynomials of degree < 3. 

3 3 2^ ■ f\ 



Remark 6 // fix) = xl then r(/) = i,M(/) = ^, and /(/) = i ^ 
iM(/) + |T(/)^/(/). 

Hence L is noi exaci in general for polynomials of degree < 4. 
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3.1 Other Points on Boundary of n Cube 
3.1.1 n = 2 

Let D = unit square. The formulas above use the vertices of D to generate 
T{f). Here we use the points (a, 0), (0, 6), (c, 1), (1, rf) on d{D), < a < 
1, < 6 < 1, < c < 1, < d < 1, to define 

T{f) = ^(/(a, 0) + /(O, b) + fie 1) + /(I, d)) 

As earlier, 

M(/) = /(l/2,l/2) 

For < A < 1, 

= AM(/) + (1 - A)r(/) 

The idea is to choose a, b, c, d and A so that Lx is exact for all polynomials 
of degree < 3. 

• La(x) = iA + |(l- A)(a + c+l) 
. L,{y) = lX + l{l-X){b + l + d) 

. Lx{x') = IX + 1{1 - X){a^ + + 1) 

. L,(y2) = iA+i(l-A)(62 + l + d2) 

. Lx{xy) = lX + l{l-X){c + d) 

. L;,(x3) = lA + i (1 - A) (a^ + c3 + 1) 

. L,(y')^lX+l{l-X)(b' + l + d') 

• L,{x'y) = lX + l{l-X){c' + d) 

• L,{xy')^lX+l{l-X){c + d') 
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Letting /(/) = f{x,y)dydx, wc have /(x) = I{y) = |, /(x^) = 

= i, I{xy) = I{x^) = I{y^) = i, I{xhj) = lixy') = i, /(x^y^) = 
I I{x^y) = I{xy^) = |, I{x^) = I{y^) = \. Setting = /(/) for the 

nine monomials above yields a system of nine polynomial equations in nine 
unknowns. Maple gives the Grobner Basis 

{3A + 6ci - 2 - 6Ad - Qd? + 6Aci^ a- d,-l + h + d,-l + c + d} 

a — d, b — 1 — d, c — 1 — d, X = f^irffj^- We are free to choose d to 
force exactness for additional monomials. However, it is not possible to get 
exactness for all fourth degree polynomials. Say we want Lx{x^y) = I{x^y) = 
\ =^ ~^ ~^S-M+t^'^' = I ^ d = 0,d = l,d = I. d = or d = l uses the 
vertices of D-the case n = 2 of the general method discussed earlier for the 
n cube. For d — ^ we have A = |. with this choice, is also exact for 
f{x,y) = xy^. Lx{x^) = ^ 7^ -^(^^) — so we do not get exactness for all 
polynomials of degree < 4. We summarize 

CR4: L(/) = i/(l/2,l/2) + l(/(l/2,0) + /(0,l/2) + /(l/2,l) + /(l,l/2)) 

is exact for all polynomials of degree < 3. Note that L(f) is also exact for 
x^y and xy^. 

Remark 7 It is not possible to choose d so that A = 0, since 6d^ — 6d + 2 
has no real roots. 

Remark 8 For the n cube in general, n >3, the natural extension would be 
to use the center of mass of the faces of Cn to generate T{f). We leave the 
details to the reader. 

Remark 9 Since C„ = [0, 1] x ■ • ■ x [0, 1], one could also use Simpson's Rule 
as a product rule. However, this does not yield any of the rules obtained 
in this section. 



11 



3.2 Connection with Interpolation 

Let 



p{x, y) - Aix^ + A2y^ + A^x + A^y + A 



Given 0<a<l,0<6<l,0<c<l,0<d<l, and /(x, y), let 

Pi = (a, 0), P2 = (0, 6), P3 = (c, 1), P4 = (1, d). We wish to choose the Aj 
so that 

p{Pk) = f{Pk), k = l,...,5 

( a" 


where P5 = (1/2, 1/2). The corresponding coefficient matrix is B — 



If a = 1/2, h = 1/2, c = 1/2, d = 1/2, then det(5) 




62 
1 

d? 

\ 4 4 
^ 7^ 0. In this 



case there is a unique interpolant, and B 



( M \ 

A2 

^3 

A, 

\ ^5 ; 



= B- 



( h \ 

b2 
h 
64 

V ^5 ; 



/ 262+264-465 \ 

261 + 263 - 46.5 
-362 - 64 + 465 
-36i - 63 + 465 

61 + 62 - 65 





1 
4 
1 
1 
4 



16 





1 

4 
1 

1 

\ 
4 



i 1 \ 
\ 1 



1 1 

' 1 



Then 



2 ly 

where 61 = /(a,0),62 
1 fi/ 



/(0,6),63 = /(c,l),64 = /(l,rf),65 = /(1/2, 1/2). j'^^A.x^ + A^y^ + A^x + 
A^y + ^5)^x^1/ = + A2) + 1(^3 + A4) + A = 

i(262 + 264 - 465 + 26i + 263 - 465) + i(-362 - 64 + 465 - 36i - 63 + 465) + 
61 + 62 - 65 = 

I62 + I64 + §65 + |6i + \h = 1/(1/2, 1/2) + |(/(l/2, 0) + /(0, 1/2) + 
/(l/2,l) + /(l,l/2)) 

= L{f) above. Thus L{f) = \f {1/2,1/2) + |(/(l/2,0) + /(0, 1/2) + 
/(1/2, 1) + /(1, 1/2)) does arise as the integral of a unique interpolant of the 
form Aix'^ + A2y'^ + ^3^ + A^y + A^. 

If a = 1, 6 = 0, c = 0, (i = 1, however, then det(-B) = 0. Thus, in this 
case, the interpolant of the form Aix"^ + A2y'^ + A^x + A^y + A^ is not unique. 

It is also natural to ask whether the cubature rule 



a 

c 
1 

1 1 



1 \ 



6 
1 
d 



1 
1 
1 

1/ 
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Tif) = i(/(l/2, 0) + /(0, 1/2) + /(l/2, 1) +/(1, 1/2)) equals the integral 
of a unique second degree interpolant to / of the form 



p{x, y) = Aixy + A^x + + 



We wish to choose the Aj so that 



The corresponding coefficient matrix is A = 



and 



f a 1\ 
6 1 
c 1 c 1 

\d d 1 1 y 

det{A) — cab — ac — cdb — dab + dac + db. Note that iia — b — c — d 
then 

det(y4) = 0. Also, if a = 1, 6 = 0, c = 0, c? = 1 (gives the vertices), then 
again dei{A) = 0. In either case, there is not a unique interpolant, p{x, y), in 
general, and thus T{f) does not arise as the integral of a unique interpolant 
of the form Aixy + A2X + A^y + A4. 



1 

2' 



4 Polygons in the Plane 

We have aheady discussed cubature formulas over triangles as a special case 
of the n simplex, and with points on the boundary other than the vertices. We 
also discussed cubature formulas over the unit square, as a special case of the 
n cube, and with points on the boundary other than the vertices. If n > 4, 
and T(/) is generated using the vertices, then the weighted combination 
AM(/) + (1 — A)T(/) is a poor CR. We examine the special case n = 6. 



4.1 6 sided Regular Polygons 

Let D be the regular hexagon with vertices 

(1 + ^3, 0), (1, 1), (-1, 1), (-l-\/3, 0), (-1, -1), (1, -1). Then the Center 
of Mass of 

D = (0,0). Using 

fix, y)dA = !\ !\ fix, y)dydx + /^^ /(^T-^avVs^ /(^' ?/)^^^^ + 

/ri'_V3 i-i^.l\wi)%f{^^ y)dydx 
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we list some useful integrals: 
Area(L)) = J^dA = A + 2^3 

JjjxdA = JoydA = J^x^dA = Joy^dA = JoxydA = JnX^ydA = 
Jj^ xy^dA = 

Xd x^dA = f + 3a/3, y^dA = | + |a/3. Let 
M(/) = (+2V3)/(0,0), 

T(/) = (4 + 2^)|(/(1 + V3, 0) + /(1, 1) + /(-1, 1) + /(-I - v^, 0) + 
/(-1,-1) + /(1,-1)) 

Lx{f) = AM(/) + (l — A)T(/). It follows easily that Lx is exact for l,x,y, 
for any A. 

Lx{x^) = (4 + 2^3) (I - |A) (^(l + v^)' + 4 + (-1 - x/s)') , and 
L\{y'^) = 4 (^4 + 2v^^ — |A^. To make Lx exact for |/^, say, we need 
4(4 + 2V3)(i-lA) = | + |V3^A=-4±||5 = i^. But then 

L,(x2) = (4 + 2V3) (i - iA) (^{l + ^y + 4+ (-1 - ^)') = 

|-\/3 + 5 7^ /(a;^) = ^ + 3-\/3. Hence L;^ is only exact for the class of 
linear polynomials. 

5 Irregular Quadrilaterals 

Let D C he the trapezoid with vertices {(0, 0), (1, 0), (0, 1), (1, 2)}(call 
them Pj). 

Jj^ dA = |, xdA = |, Jjy ydA = | ^ Center of mass of D is (5/9, 7/9). 
Define the following linear functionals: 

M(/) = Area(D) /(5/9, 7/9) = ^/(5/9, 7/9) 

T(/) = Area(L')( J ^ /(P,)) = |(/(0, 0) + /(I, 0) + /(0, 1) + /(I, 2) 
For fixed A, < A < 1, 

Lx = AM(/) + (1 - A)r(/) 

To make Lx exact for x, we have Lx{x) = ^A + | = |^A = 1. Then 
-^a(|/) = ^-^ + 1 = 1- However, LA(a;|/) = f 7^ Jq Jq^^ xydydx = ^. So 
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using the vertices of d{D) for T{f) only gives exactness for linear functions 
in general. 



Remcirk 10 It is interesting to see if there are any weights wi, W2, w^, W4, W5 
such that the CR L{f) = ^1/(5/9, 7/9) + ^2/(0, 0) + ^3/(1, 0) + ^4/(0, 1) + 
w^f (1,2) is exact for all polynomials of degree < 2. Above we considered the 



special case where the weights would have the form wi = ^ X, w 



A), 



2, 3, 4, 5. Here is the resulting linear system: wi + W2 + ws + W4^ + W5 ~ 



Wi + Ws + W5 = \, \wi + W4 + 2W5 



7 2'i 



J_ 

12' 



2' 



^wi + 'u;4 + 4^5 = f , ^wi + 2^5 = ^ 



81 



The augmented matrix is A 



( 1 

5 



9 

25 

81 
49 

§1 
\ 35 
\ 81 



24 
1 









1 

1 



1 






1 



1 



1 





2 

5 

f 

6 

J_ 

12 
5 

iz 

24 



and the RREF of 



24 / 



%S 



( 1 








V 






1 












1 






\ 







1 ) 



Hence the system has no solution, and thus there 



are no weights wi,W2,W3,W4,W5 such that the CR L(f) — wi/(5/9, 7/9) + 
W2f (0,0) + wsf (1,0) + w J{0,1) + W5f (1,2) is exact for all 

polynomials of degree < 2. One can, however, choose weights so that 
L{f) is exact for all polynomials of degree < 2, except for xy, say. The 
corresponding augmented matrix is 



( 



81 
49 

81 



1 








12 

5 



which has row echelon form: 



( 1 










1 











1 











1 











1 



120 

29 

240 
31 

240 
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5.0.1 Other points on Boundary of Trapezoid 

We will try using other points on d{D) to generate T(/)— say (a, 0), (0, b), (1, c), {d, d+ 
1), with 0<a<l,0<6<l,0<c<2,0<d<l. Then 

Lxif) = A^/(5/9, 7/9) + (1 - A)^(/(a, 0) + /(I, c) + /(0, 6) + /(d, d + 1) 

• LaIx) = |A + I (1 - A) (a + 1 + d) 
. Lxiy)^lX + l{l-X){c + b + d+l) 
. L;,(x?/) = ||A + I (1 - A) (c + d (rf + 1)) 
. L,(x2) = |A + I (1 - A) (a^ + 1 + d^) 

. L,iy') = ||A + I (1 - A) (c^ + 6^ + (d + 1)^) 

Wc set each of the above expressions in A, a, b, c, d equal to the corre- 
sponding integral of / over D to yield the system: 

|A + |(l-A)(a + l + d) = |,|A+|(l-A)(c + 6 + d+l) = I 
|a + |(1-A) (c + rf(ci + l)) = l|,|A + |(l-A) (a^ + l + d2) = X 

||A + |(1-A) {c^ + b' + {d+lf) = l. 

We used Maple to find the following Grobner basis: 

{392A-163, 9a+9d-ll, 816-99d+20, 180d-191+81c, -22671d+6583+18549d^} 



The last equation has solutions d = ^ ± ^v3893 ~ . 74734, . 47488 
Wc shall use d = ^ + -^^/SSdS. Solving the other equations yields 
180d - 191 + 81c = ^ c = 1 - 25^1 V3893 ^ . 69726 
816 - 99d + 20 = =^ 6 = i + 41^2 \/3893 w . 6665 
9a + 9d-ll = 0^a=i| - ^V^S93 ^ . 47488 
392A -163 = 0^A = g|Ri. 41582. We summarize 

CR5: Let D C R"^ be th e trapezoid with vertices {(0, ), (1, 0), (0, 1), (1, 2)}. 

Let a = i| - -^VS893 ^ . 47488, 6 = | + ^V^893 ^ . 6665, 

c = 1 - 2561^/3^ ^ ■ 69726, d = ^ + 4^8^1893 ^ . 74734, and A = |f ^ 
. 41582. Then L{f) = A|/(5/9, 7/9) + (1 - A)|(/(a, 0) + /(I, c) + /(o', b) + 
f{d, d + 1)) is exact for all polys, of degree < 2. 

Note that L{x^) = f||§ .44092, while I{x^) = | = .45. So L is not 
exact for cubics. 
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6 Unit Disc 



Let D be the unit disc. The n roots of unity Zk = e^'^*'^/"', k = 1, n, can be 
used to generate T(/) = ^ Zlfc=i /(cos( 27ik/n, sin(27rA;/n)). Since the center 
of mass oi D — (0, 0) and area(£)) = tt, the analogy of our formulas for 
the triangle and the square is M(/) = 7r/(0,0), Lx{f) = A7r/(0,0) + (1 - 
I]fe=i /(cos( 27rA;/n, sin(27rA;/n)). It is not hard to show, however, that 
the only good choice is n = 4. This gives the formulas 

Lx(f) = A7r/(0, 0) + (1 - A)^(/(l, 0) + /(0, 1) + /(-1, 0) + /(O, -1)) 

We list the following integrals without proof: 

• If m and n are even whole numbers, then 

/■ /■ x-v-dA = ^ ( jn-mm-iy. \ 

J Jd ^ m + n + 2 V2"»+"-3(|-l)!(f -l)!(Hi±!i)!y 

• If m is an even whole number and n = 0, then 



Id ^ m + 2 V2"^-2(f - l)!(f)!^ 

• If n is an even whole number and m — 0, then 



J x^'y^'dA = 



TT / (n-1)! 



n + 2 V2"-2(n_ 1)1(1)!^ 
If m and/or n is an odd whole number, then 



dA^O 



It is then easy to prove that the choice A = | gives 
CR6: L(/) = |/(0, 0) + 0) + /(0, 1) + /(-1, 0) + /(O, -1)) 

is exact for all polynomials of degree < 3. 
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7 Summary of Cubature Rules 



These are the CRs we derived as a generahzation of Simpson's Rule for 
functions of one variable. If D„ is a polygonal region in i?", let Po,...,Pn 
denote the n + 1 vertices of Dn and Pn+i — Center of Mass of Dn. The 
generalization is based on the weighted combination L\ = \M{f) + (1 — 
A)T(/), where M(/) = Vol(D„) /(P„+i), T(/) = yo/(D„)(;;^ Ef=o /(^.O)- 
Unless noted otherwise, all of the rules have the following property: 
All of the weights are positive, all of the knots he inside the region, and 
all but one of the knots hes on the boundary of the region. It is desirable 
to have as many knots as possible on d{Dn) if one subdivides the region and 
compounds the CR. 

7.1 n Simplex T„ 

• Using the vertices, L{f) = (^^^/(Pn+i) + j^E]=o f{Pj) is exact 
for all polynomials of degree < 2. 

• Using points other than the vertices: Letting {Qk} denote the center of 
mass of the faces of T„, L(/) = - (n - 2) n±ii./(i/(n + 1), l/(n + 

^)) + (^"2)! Sfcii fiQk) is exact for all polynomials of degree < 2. All 
but one weight is positive if n > 2. 

7.2 Unit n Cube Cn 

• General n: Let Pq, Pm-i denote the m vertices of Cn, m ~ 2"'. Then 

L(/) = f /(1/2, 1/2) + 1^ E?lo' fiPj) is exact for aU polynomials 
of degree < 3. 

. n = 2: L{f) = i/(l/2, 1/2) + l(/(l/2, 0) + /(0, 1/2) + /(1/2, 1) + 
/(1, 1/2)) is exact for all polynomials of degree < 3. 

7.3 Unit Disc D 

L{f) = f /(O, 0) + f (/(I, 0) + /(0, 1) + /(-I, 0) + /(O, -1)) is exact for aU 
polynomials of degree < 3. 



18 



7.4 A Trapezoid in the Plane 



Let D C Ei^he the tr apezo id with vertices {(0, 0), (1,0), (0, 1), (1, 2)}. 
Let a = - -^V^Q^ . 47488, 6 = | + -^V^93 ^ . 6665, 
c = 1 - ^./M93 ^ . 69726, d = ^ + ^V^^ ^ . 74734, and A = if ^ 

. 41582. Then = A|/(5/9, 7/9) + (1 - A)|(/(a, 0) + /(I, c) + /(O, b) + 

f{d, d + 1)) is exact for all polys, of degree < 2. 
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